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We consider the dynamics of systems of self propelling particles with nonholonomic constraints. 
A continuum model for a discrete algorithm used in works by T. Vicsek et al. is proposed. For a 
case of planar geometry the finite flocking behavior is obtained. The circulation of the velocity field 
is found not to be conserved as a consequence of the nonholomicity. The stability of ordered motion 
with respect to noise is discussed. An analogy with the kinetics of charges in superconductors is 
noted. 
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The emergence of ordered structures in dynamic sys- 
tems is a long standing problem in physics. Generally 
speaking one deals with dynamic phase transitions gov- 
erned either by external or internal noise. Most inter- 
esting is the arising of ordered motion caused by the in- 
ternal dynamics of the system. Recently, there has been 
a growing interest in studying the collective behavior in 
systems of self-propelling particles (SPP). One may dis- 
tinguish systems of two types: 

The first one is for systems of "unintelligent" parti- 
cles interacting via real physical forces produced by the 
background in which they are moving (micelles, bacteria, 
etc). Here the driving forces are caused by the gradients 
of chemical factors (concentrations, chem. potential etc.) 
or physical factors (light, potential and dissipative fields 
etc.), which influence the motion. In particular, the earth 
magnetic field is vital for the orientation during long dis- 
tance migration of biological species like birds or turtles 
0, ■ It is clear that the absence of conservation of 
translational and angular momentum is a direct conse- 
quence of these external factors |3| . 

The second class is formed by systems of particles, 
which interact via nonholonomic constraints imposed on 
their velocities. We would like to stress that the non- 
holonomicity of such systems clearly expresses the "in- 
telligent" nature of the particles, since such constraints 
need instant exchange of information (visual or any other 
sensorial) between the particles and their environment. 
This explains the coherent motion and arising ordered 
patterns in the dynamics of systems like crowds, traffic 
or flocking. Usual potential gradients or other physical 
forces are not relevant for the collective behavior, though 
the particles need some physical source of energy to sus- 
tain the constraints. The physical origin of the nonholo- 
nomicity is the force which acts on the particle due to its 
interaction with the background (earth, air, liquid sub- 
strate etc.). Since we are not interested in the dynamic 
degrees of freedom of the background we lose this infor- 



mation. Such loss of dynamic information leads to the 
breaking of the conservation of the (angular) momen- 
tum and effectively to the nonholonomic constraints. In 
general it means that the system with such constraints 
is not closed and therefore its dynamics is not Hamilto- 
nian, though it does not mean that the energy dissipates. 
Rather a redistribution among the dynamic degrees of 
freedom takes place. But the very form of the constraint 
is determined by the "intellect " of a particle which uses 
the information about the environment and moves ac- 
cordingly. Note that the numerical algorithm used in 
Ref. |j] modelled that kind of systems. For shortness we 
call it the Czirok- Vicsek automaton or algorithm (CVA). 
There are also modifications of this algorithm, which dif- 
fer from the CVA by inclusion of potential interparticle 
forces 0, external regular and stochastic fields p. We 
will consider CVA as the minimal model for the collective 
behavior since the main cause for the self-organization in 
the system is the nonholonomicity of CVA dynamic rule. 

To our knowledge none of the hydrodynamic models 
proposed earlier is aimed to reflect the essence of the 
CVA - its conservative and distinctly nonholonimic char- 
acter. To a great extent such models are modifications 
of the Navier-Stokes equation. Such an approach is cer- 
tainly valid for microorganisms floating in a medium. It 
is, however, hardly adequate for the collective behavior 
of "intelligent" boids, e.g. birds or drivers in traffic flow 
7], when the individual behavior is determined rather 
by the instant exchange of information with the envi- 
ronment than by the action of some interparticle forces. 
The key point here is the nonholonomicity of the system. 
In addition all the models based on the modifications of 
the Navier-Stokes equation like Ref. 0, E| include addi- 
tional phenomenological terms which generate the "spon- 
taneous" transition to the state of ordered motion. Such 
terms are added much in formal analogy with Landau's 
theory of equilibrium continuous phase transitions and 
the notion of an order parameter but any firm base for 
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such an analogy is not given. Thus the ordered state in- 
troduced in such a model is rather an ad hoc assumption 
than the natural consequence of the underlying interpar- 
ticle interactions. The interpretation of the viscous term 
for such systems is also completely unclear. 

In Ref. Q it was noted that the CVA can be considered 
as the dynamic XY-mode\. It was concluded also that 
the cause of the ordering is the convective term. The 
AF-model is a Hamiltonian dynamic system. The latter 
property is the key point since for Hamiltonian dynamic 
flow there is an ergodic measure. Due to nonholonomic 
character the dynamic rule of the CVA even in the static 
limit |v| — > 0, i.e. when there is no any transfer flow 
term like v • V, has no canonical Hamiltonian form. The 
general fact is that such systems do not have an invari- 
ant measure with respect to the dynamic flow. In fact 
the ordered state in the CVA appears as the fixed point 
(attractor) of its averaging dynamic rule. Note that from 
the point of view of the theory of dynamic systems the 
compactness of the phase space (i.e. the space of posi- 
tions and velocities) is very important so the question 
about influence of the boundary conditions usually used 
in simulations on the ordering the system also needs dis- 
cussion. 

From such a nonholonomic point of view the appear- 
ance of ordered motion in the CVA and similar dynamic 
systems is a trivial consequence of their nonholonomic- 
ity. Indeed, the breaking of conservation of the (angular) 
momentum is due to nonholonomic constraints which as 
has been mentioned above mean that the system is not 
closed. 

Here we discuss a hydrodynamic model which can be 
considered as the continuous analogue of the discrete dy- 
namic automaton proposed in Ref. |J] for the SPP sys- 
tem. It manifestly takes into account all local conser- 
vation laws for the number of particles and the kinetic 
energy. The self propelling force and the frictional force 
are assumed to balance each other. 

The algorithm used in Ref. corresponds to the fol- 
lowing equation of motion of a particle: 

-^V; =U f X Vj, (1) 

where uii is the "angular velocity" of i-th boid, which 
depends on what happens in the neighborhood. It is 
assumed also that the number of particles is conserved. 
Like algorithm in Ref. , Eq. Q distinctly expresses the 
conservation of the kinetic energy. 

From a physical point of view it is natural that the 
hydrodynamic model corresponding to the CVA is based 
on the following equations: 



± J n(r,t)dV = 0, 



v 
d 
dt 



jiv 2 dV = , 



(2) 



(3) 



where n(r,t) and v(r, t) are the number density and the 
Eulerian velocity. The volume V moves along with the 
velocity field. The first condition is the conservation of 
number of particles. As usual we can rewrite this condi- 
tion in the differential form: 



dn 
~dt 



div (nv) = 0. 



(4) 



The second constraint Eq. © means that the kinetic 
energy of a Lagrange particle is conserved, i.e. 



4- fnv 2 dV 
dtj 

v 



v" I — h ndivv ! dV 



n——dV = 0. 

dt 



(5) 



The first integrand vanishes due to the conservation of 
the particle number. As a consequence the second inte- 
gral also vanishes for an arbitrary choice of the volume 
V which in view of the natural condition n > leads to: 



dt 



|v(r,t)| 2 =0. 



(6) 



This implies that a pseudovector field u> exists such that 



dt 



v (r, t) = u (r, t) x v (r, t) . 



(7) 



This equation, which is the continuous analogue of 
Eq. 0}, nas now been derived from the conservation of 
particle number and kinetic energy, Eq. and Eq. ©. 
We will model u) by: 

u>(r,t) = J K (r - r') n(r', t) rot v(r', t) dr' , (8) 

which has the proper pseudovector character and heuris- 
tically may be considered as the continual analog of the 
CVA dynamic rule. There are other possible choices like: 



u>(r,t) = J K(r-r') Vn(r',i) x v(r',i)dr' 



(9) 



and combinations of the two. The averaging kernels K 
and K should naturally decrease with the distance in 
realistic models. We concentrate our discussion on the 
case Eq. |JSJ. Equations (0J and (7J obviously have the 
uniform ordered motion with a constant both density and 
velocity as trivial solution. One may consider Eq. Q as 
the equation of motion of a charge in a magnetic field, 
where v is the charge velocity and u> is proportional to 
the magnetic field. One may even include an "electric 
field" and dissipative (collision) terms in Eq. @: 



dt 



v(r,i)= f + w(r,i)xv(r,t)-£v(r,i) . (10) 
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where £ _1 is the mean free time. We will not consider 
such an extended model here since algorithm used in |4| 
is distinctly conservative. Further we exploit the anal- 
ogy of Eq. O and the equation of motion for charges 
in superconductors (see e.g. Ref. 0). In our model a 
current density corresponding to the particle velocity at 
some point r depends on "magnetic field" u>, i.e. "vec- 
tor potential", at all neighboring points r' within some 
region of coherence. Therefore the relations © or 
can be considered as a corresponding nonlocal relations 
between the current density and the vector potential in 
nonlocal Pippard's electrodynamics of the superconduc- 
tors. The situation in the system under consideration 
is more complex. In electrodynamics of superconductors 
the external magnetic field is the main cause of the vor- 
tical motion of charges since their own magnetic field is 
negligibly small and does not lead to formation of the 
ordered motion. In our case the " magnetic field" u) itself 
depends on the motion of "charges", i.e. particles, and 
vise versa, which leads to a nonlinearity of the system. 
Depending on the parameters one can expect either the 
direct current state since the system is conservative or 
vortical states like Meissner currents or Abrikosov vor- 
tices. 

Using such an analogy let us find the conditions for the 
existence of the stationary vortical states. 
We can rewrite Eq. 10 in the following form: 



— + — = rot(vxW) 



(11) 



where W(r, t) = rot v — ui . Thus it follows that if W(r, t) 
is equal to 0, then duj/dt = and therefore u> = rotv 
is independent of the time. Such states are naturally in- 
terpreted as stationary translational u) — or rotational 
u> ^ regimes of motion. For the model © together 
with W = we get the integral equation: 



J K (r — r') n(r') rot v(r') dr' = rot v(r) , 



(12) 



which determines such states. Equation 1)12(1 gives sta- 
tionary vortical motion, represented by the vector field 
with |v| = const. From here it follows that the vortic- 
ity of the velocity field is an eigenstate of the integral 
operator with n(r) as the corresponding weight factor. 
It should be noted that these stationary states do not 
exhaust all stationary vortical states, since in general 
Vv 2 ^ 0. 

We further scale K by multiplying with some n* and 
similarly scale the density by dividing by n* . Further- 
more we restrict our discussion to the simple case of a pla- 
nar geometry with averaging kernel in (JSJ) as ^-functional: 

K(r - r') = s 5(r - r') , s = ±1. (13) 
We will call this the local model. For this case Eqs. (0J 



and 0) take the form: 
On 



— + div (nv) = 0, 
d 

—v = s n rot v x v . 

dt 



(14) 
(15) 



These equation <|15[l can be obtained as a special case 
of the corresponding one in Ref. 8]. In our work the 
foundation of the terms is given. We do not use general 
symmetry arguments to take all terms of a certain sym- 
metry into account irrespective of their physical mean- 
ing. The corresponding local model with (5-kernel for © 
may be identified with rotor chemotaxis (i.e. caused by 
chemical field) force introduced in Ref. [llj if one takes 
into account simple linear relation between field of food 
concentration and the number density of boids, which is 
surely valid for low concentrations of food and bacteria. 

The parameter s of the local model l|15fl distinguishes 
different physical situations concerning the microscopic 
nonholonomical constraint. To see this we find the sta- 
tionary radially symmetric solutions of ((14JI and H15J) ■ As 
usual, we search for the solutions of the form n = n(r), 
v = v v (r)e,p. The continuity equation l|14l) is satisfied 
trivially. Substituting this into Eq. l(TK|) we finally ob- 
tain: 



/ \ Cst 



1 



r n(r) 



dr 



(16) 



where ro is the cut-off radius of the vortex. This is the 
core of the vortex. The constant C s t is determined by 
the circulation of the core § v d\ = C s t ■ 

r=r 

The spatial character of the solution strongly depends 
on the parameter s. If s = — 1 the infinitely extended 
distributions for n(r) are allowed, e.g. n(r) oc r~ a , a > 

0. They lead to localized vortices with exponential decay 
of angular velocity. If s = +1 only compact distributions, 

1. e. n(r) = outside some compact region, are consistent 
with the finiteness of the total kinetic energy i.e. they 
corresponds to finite number of particles J n dV < oo . 
As an example we may give: 



n(r) = 




r <r < R, 
otherwise 



(17) 



Substituting Eq. (JHI) into Eq. JTSJ) one obtains 

Cst 



27rr 



exp 



ro 



1 



arctanh 

R 



ro- 1 

(18) 

The corresponding component of the velocity v = 
for such a case is shown on Fig. ^ 

Taking the rotation on both sides of Eq. I|15|) one ob- 
tains the following equation for the vorticity in a case of 
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FIG. 1: Velocity V*(r/R) = 2%Rv^{r)/Cst in the vortex of 
the local model with n(r) given by Eq. JT7) at r /R = 1/3. 



planar geometry: 

dv f dn 

rot — = — s — — rot v — n (v • V) rot v 



(19) 



This implies that for s = +1 the vorticity is damped 
by compression along the flow and therefore such a flow 
is stable with respect to the vortical perturbations. For 
s = — 1 the vorticity is damped by expansion. 

The hydrodynamics of nonholonomic fluid under con- 
sideration differs essentially from the potential dynamics 
of ideal fluids |l2( due to the nonpotential character of 
the equation of motion (JJJ . In ideal fluids the nontrivial 
stationary flows are only possible under the influence of 
an external force (pressure gradient or reaction of bound- 
aries). Moreover the vorticity in the nonholonomic fluid 
is not coserved. This is not a case for usual ideal fluids 
where the vorticity is the integral of motion. 

The first term in right-hand side of Eq. i|19|) shows the 
influence of compression on the evolution of the vorticity. 
The second term represent the modified spatial transfer 
of the vorticity along the flow. Since this term can be ex- 
cluded locally by the choice of the instant local reference 
frame which moves accordingly with the flow we consider 
the first term as the main source of the vorticity. In such 
an approximation we can write: 



dv dn 
rot — — = — s— - rot v. 

dt dt 



The circulation is defined by 



C(t) = <t> vdl = / rotvdS, 



(20) 



(21) 



where the integration contour and the corresponding sur- 
face area move along with the velocity field. The time 
derivative of the circulation is: 



dt 



C 



dv 

rot — dS . 

dt 



(22) 



Thus the circulation does not conserve in contrast to the 
ideal fluid model, where the microscopic interactions are 
of holonomic character. 

The total momentum P = J n v dV does not conserve 
too. For the local model we can write: 



dt 



' rot v x v dV . 



(23) 



From Eq. (|23[1 it follows that the damping of vortical 
part of the velocity leads to the formation of the state of 
uniform motion with P = const. 

Here we consider the influence of noise on the stability 
of the flow with respect to the vortical perturbation. It 
is clear that instability with respect to such perturba- 
tions drives the system to disordered state. Inclusion of 
stochastic noise can be done in a way analogous to that 
used in Ref. Q: lj (r, t) = uj (r,t) + 5u)(r,t) , where 
loq = sn rot v is the same contribution as before and 5uj 
is the stochastic contribution. These fluctuations lead 
to fluctuations of the density and velocity fields. Re- 
placing dn/dt by an average value 1/r plus a fluctuating 
contribution SL(t) in Eq. (|20() one obtains for the above 
described local model with s = +1: 



dt 



i + SL ) C . 



(24) 



The simplest model for the noise is the Gaussian white 
noise approximation: 

(6L(t)6L{t')) = 2TS(t-t'). (25) 

The stochastic equation l|24|) has the solution: 

t" 



C{t) = C exp 



exp(-W(t)) , 



(26) 



where W(t) = jSL(t')dt' is the Wiener process |13j . 
o 

Averaging over the realization of the stochastic process 
we get the averaged evolution of the vorticity: 



.on, =r„^q»( - J 



i-rr' 



(27) 



where f is the relaxation time of the circulation in the 
system. For large enough noise r T > 1 the system be- 
comes unstable. For r T < 1 the system is stable and the 
circulation decays to zero. When tT increases to 1 the re- 
laxation time f goes to infinity, a result similar to critical 
slowing down near the critical point. These estimates are 
modified when other terms, which have been neglected, 
are taken into account but we believe that qualitatively 
the obtained results remain unchanged. 

It should be noted that the average kernels can also 
contain noise contributions. In view of the above results 
for the character of stationary states and vorticity relax- 
ation, which depends on the sign of parameter s, this 
case needs a more thorough investigation. 
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In conclusion we have constructed a continuum SPP 
model with particle number and kinetic energy conser- 
vation. We found in 2D that vortical solutions exist 
for the model and that they show a finite flocking be- 
havior. These solutions, which qualitatively reproduce 
some results of Ref. Q , were obtained without imposing 
any boundary conditions on the velocity field. The non- 
holonomic constrains were found to lead to a circulation 
which was not conserved. The influence of noise on the 
stability of the system was discussed. 
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